function [difff, varargout] = Steady(rr_guess,varargin)  

Par = varargin{1};
Step  = varargin{2};
Trans = varargin{3};
QA_init = varargin{4};

w_adj = Par.lmd^(-Par.mbar/(1-Par.bet))*((1+Par.kapa)*(1+Par.trf_A))^(-1/Par.bet)+10^-3; % Ensure m_M < m_bar: Not everything imported!

rr_guess    = exp(rr_guess);

qrat_A = rr_guess(3);
qrat_B = rr_guess(4);

Qrat = qrat_A/qrat_B;

if Qrat<w_adj; qrat_B = qrat_A/w_adj; Qrat = w_adj; end        
                
wq_A = (1-Par.bet)*Par.eta^(Par.bet-1)*qrat_A^(-Par.bet);
wq_B = (1-Par.bet)*Par.eta^(Par.bet-1)*qrat_B^(-Par.bet);    

xk_A = (Par.pi_frn/Par.bet)^(1/(1-Par.bet))/wq_A^(1/Par.bet)*Par.L_B;
xk_B = (Par.pi_frn_toUS/Par.bet)^(1/(1-Par.bet))/wq_B^(1/Par.bet)*Par.L_A;

Par.mreal_X = (1-Par.bet)*(+1/Par.bet*log((1+Par.kapa)*(1+Par.trf_B))-log(Qrat))/log(Par.lmd); % possibly positive
Par.mreal_M = (1-Par.bet)*(-1/Par.bet*log((1+Par.kapa)*(1+Par.trf_A))-log(Qrat))/log(Par.lmd); % possibly negative

Par.madj_X = round(Par.mreal_X)-1; % largest step with zero exports in A
Par.madj_X = min(max(Par.madj_X,-(Par.mbar+1)),Par.mbar);

Par.madj_M = -(round(Par.mreal_M)+1); % smallest step with zero imports in A (absolute val)
Par.madj_M = min(max(Par.madj_M,-(Par.mbar+1)),Par.mbar);

Par.mcut_X = Par.madj_X+1;
Par.mcut_M = Par.madj_M+1;

eps_X = 1-(Par.mreal_X-round(Par.mreal_X)+0.5);
eps_M = -(Par.mreal_M-round(Par.mreal_M)-0.5);

madj_XB = round((1-Par.bet)*(1/Par.bet*log((1+Par.kapa)*(1+Par.trf_A))+log(Qrat))/log(Par.lmd))-1;
madj_XB = min(max(madj_XB,-(Par.mbar+1)),Par.mbar);

Step{5} = [];

Step{5}(:,1) = [ones(Par.mbar-Par.mcut_X,1)*(Par.pi_dom*Par.L_A+Par.pi_frn*Par.L_B);
                ones((Par.mbar-Par.madj_X)>0,1)'*Par.pi_dom*Par.L_A+eps_X*Par.pi_frn*Par.L_B;
              ones(Par.madj_X+Par.madj_M+1,1)*Par.pi_dom*Par.L_A; 
              ones((Par.mbar-Par.madj_M)>0,1)'*eps_M*Par.pi_dom*Par.L_A;
              zeros(Par.mbar-Par.mcut_M,1)]*wq_A^(1-1/Par.bet);
Step{5}(:,2) = [ones(Par.mbar-Par.mcut_M,1)*(Par.pi_dom*Par.L_B+Par.pi_frn_toUS*Par.L_A);
                ones((Par.mbar-Par.madj_M)>0,1)'*Par.pi_dom*Par.L_B+(1-eps_M)*Par.pi_frn_toUS*Par.L_A;
              ones(Par.madj_X+Par.madj_M+1,1)*Par.pi_dom*Par.L_B; 
              ones((Par.mbar-Par.madj_X)>0,1)'*(1-eps_X)*Par.pi_dom*Par.L_B;
              zeros(Par.mbar-Par.mcut_X,1)]*wq_B^(1-1/Par.bet);
          
Step{6}(1,:) = [ones(1,Par.mbar+1+Par.madj_M) ones((Par.mbar-Par.madj_M)>0,1)*eps_M zeros(1,Par.mbar-Par.mcut_M)]; 

Step{6}(2,:) = [ones(1,Par.mbar+1+Par.madj_X) ones((Par.mbar-Par.madj_X)>0,1)*(1-eps_X) zeros(1,Par.mbar-Par.mcut_X)]; 

steps      = Step{1};
Trans_innov = Step{2};
Trans_fail = Step{3};
Trans_exo  = Step{4};
prft_fnc_A = Step{5}(:,1);
prft_fnc_B = Step{5}(:,2);
wage_dom_A = Step{6}(1,:);
wage_dom_B = Step{6}(2,:);

Param = structvars(Par);
for ctr0 = 1:size(Param,1); eval(Param(ctr0,:)); end

Transit = structvars(Trans);
for ctr0 = 1:size(Transit,1); eval(Transit(ctr0,:)); end

% ===== Steady State Variables =====

vv_init = log([prft_fnc_A(1)/rr_guess(1)*ones(mtot,1);
               prft_fnc_B(1)/rr_guess(2)*ones(mtot,1)]);
           
% if exist('vv_guess.mat','file')>0; load vv_guess; vv_init = log(vv); end

diff_check = 0; 
if diff_check
    
    load rr_diff;
    assist = max(abs(difff))<5*10^-2;
    
else 
    
    assist = 0;
    
end

if ~assist
    
    [vv_ss,exit_ss] = Steady_set_iter(vv_init,Par,Step,rr_guess);
    vv = vv_ss;
    vv_init = log(vv_ss);
    
end

if assist || exit_ss<1
    
    if assist 
    
        display(' ')
        display(['          Getting assistance for SS !!!']);
        display(' ')
        
    else        
    
        display(' ')
        display(['          Getting assistance for SS - ' ...
                'NO convergence with  Steady_set_iter !!!']);
        display(' ')
        
    end
        
    opt_steady = optimset('Display','final');
    
    steady_new = @(N)Steady_set(N,Par,Step,rr_guess);
    
    [vv_ss,fval_ss,exit_ss,output_ss] = fsolve(steady_new,vv_init,opt_steady);
    
    if exit_ss < 1 && exist('vv_guess.mat','file')>0;
        
        display(' ')
        display(['          Getting assistance for SS - another guess !!!']);
        display(' ')
        
        load vv_guess
        
        vv_init = log(vv);
        
        opt_steady = optimset('Display','final','MaxFunEvals',17600,'MaxIter',1600);
        
        [vv_ss,fval_ss,exit_ss,output_ss] = fsolve(steady_new,vv_init,opt_steady);
        
    end

    vv = exp(vv_ss);

end

% if exit_ss==1; save vv_guess vv; end


% ===== Value Functions  and R&D decisions =====
    
vv_iA = vv(1:length(vv)/2);
vv_iB = vv(1*length(vv)/2+1:2*length(vv)/2);

xx_iA = ((Trans_innov'*vv_iA-...
        vv_iA)/(alfa_A*(1-tau_A))).^(1/(gama_A-1));
xx_iB = ((Trans_innov'*vv_iB-...
        vv_iB)/(alfa_B*(1-tau_B))).^(1/(gama_B-1));

xx_eA = ((Trans_innov'*vv_iA...
            -0)/alfa_eA).^(1/(gama_eA-1));
xx_eB = ((Trans_innov'*vv_iB...
            -0)/alfa_eB).^(1/(gama_eB-1));
        
% ===== Transition Matrices =====

Qtr_A_idom = Qloc_idom*diag(xx_iA);
Qtr_A_ifrn = Qloc_ifrn_L*(diag(xx_iB)*Qloc_frn_R);
Qtr_A_edom = Qloc_edom*diag(xx_eA);
Qtr_A_efrn = Qloc_efrn_L*(diag(xx_eB)*Qloc_frn_R);

Qtr_B_idom = Qloc_idom*diag(xx_iB);
Qtr_B_ifrn = Qloc_ifrn_L*(diag(xx_iA)*Qloc_frn_R);
Qtr_B_edom = Qloc_edom*diag(xx_eB);
Qtr_B_efrn = Qloc_efrn_L*(diag(xx_eA)*Qloc_frn_R);

Trans_matA = (Qtr_A_idom+Qtr_A_ifrn+Qtr_A_edom+Qtr_A_efrn+Qtr_exo)*dt;
Trans_matB = (Qtr_B_idom+Qtr_B_ifrn+Qtr_B_edom+Qtr_B_efrn+Qtr_exo)*dt;
          
QB_init = diag(lmd.^steps)*flip(QA_init);

QA_next = 0; QB_next = 0;
          
for ctr1 = 1:time_ss
    
    QA_init = QA_init*(ctr1==1)+QA_next*(ctr1>1);
    QB_init = QB_init*(ctr1==1)+QB_next*(ctr1>1);
    
    QA_next = QA_init+Trans_matA*QA_init;
    QB_next = QB_init+Trans_matB*QB_init;
    
    IA_init = prft_fnc_A'*QA_init-...
            alfa_A/gama_A*xx_iA'.^gama_A*QA_init-...
            alfa_eA/gama_eA*xx_eA'.^gama_eA*QA_init+...
            pi_dom/(1-bet)*wq_A^(1-1/bet)*L_A*wage_dom_A*QA_init+...
            pi_frn_toUS/(1-bet)*wq_B^(1-1/bet)*L_A*flip(1-wage_dom_A)*QB_init+...
            wq_A*sum(QA_init)+Par.rev_US*trf_A*((1+trf_A)*(1+kapa)*eta/(1-bet))*wq_B*xk_B*flip(1-wage_dom_A)*QB_init;  
    IB_init = prft_fnc_B'*QB_init-...
            alfa_B/gama_B*xx_iB'.^gama_B*QB_init-...
            alfa_eB/gama_eB*xx_eB'.^gama_eB*QB_init+... 
            pi_dom/(1-bet)*wq_B^(1-1/bet)*L_B*wage_dom_B*QB_init+...
            pi_frn/(1-bet)*wq_A^(1-1/bet)*L_B*flip(1-wage_dom_B)*QA_init+...
            wq_B*sum(QB_init)+Par.rev_FN*trf_B*((1+trf_B)*(1+kapa)*eta/(1-bet))*wq_A*xk_A*flip(1-wage_dom_B)*QA_init;     

    IA_next = prft_fnc_A'*QA_next-...
            alfa_A/gama_A*xx_iA'.^gama_A*QA_next-...
            alfa_eA/gama_eA*xx_eA'.^gama_eA*QA_next+...
            pi_dom/(1-bet)*wq_A^(1-1/bet)*L_A*wage_dom_A*QA_next+...
            pi_frn_toUS/(1-bet)*wq_B^(1-1/bet)*L_A*flip(1-wage_dom_A)*QB_next+...
            wq_A*sum(QA_next)+Par.rev_US*trf_A*((1+trf_A)*(1+kapa)*eta/(1-bet))*wq_B*xk_B*flip(1-wage_dom_A)*QB_next;  
    IB_next = prft_fnc_B'*QB_next-...
            alfa_B/gama_B*xx_iB'.^gama_B*QB_next-...
            alfa_eB/gama_eB*xx_eB'.^gama_eB*QB_next+...
            pi_dom/(1-bet)*wq_B^(1-1/bet)*L_B*wage_dom_B*QB_next+...
            pi_frn/(1-bet)*wq_A^(1-1/bet)*L_B*flip(1-wage_dom_B)*QA_next+...
            wq_B*sum(QB_next)+Par.rev_FN*trf_B*((1+trf_B)*(1+kapa)*eta/(1-bet))*wq_A*xk_A*flip(1-wage_dom_B)*QA_next; 
    
    QwA_next = wage_dom_A*QA_next+flip(1-wage_dom_B)*QA_next*(1+kapa)/((1+kapa)*(1+trf_B))^(1/bet);
    QwB_next = wage_dom_B*QB_next+flip(1-wage_dom_A)*QB_next*(1+kapa)/((1+kapa)*(1+trf_A))^(1/bet); 
    
    gA_ss_vec(ctr1) = (IA_next-IA_init)/IA_init;
    gB_ss_vec(ctr1) = (IB_next-IB_init)/IB_init;
    
    QQbar_rat(:,ctr1) = [sum(QA_next)/QwA_next; sum(QB_next)/QwB_next];
    
end

% ---- Check if SS reached ------

Tol_g = 10^-3;
iter1 = 1;
g_rat = 1;

g_check = max(abs(((1+gA_ss_vec(end-4/dt))^(1/dt)-1)/...
            ((1+gA_ss_vec(end))^(1/dt)-1)-1),...
            abs(((1+gB_ss_vec(end-4/dt))^(1/dt)-1)/...
            ((1+gB_ss_vec(end))^(1/dt)-1)-1))>Tol_g;

while g_check && g_rat>Tol_g
    
    QA_init = QA_next;
    QB_init = QB_next;
    
    QA_next = QA_init+Trans_matA*QA_init;
    QB_next = QB_init+Trans_matB*QB_init;
        
    IA_init = prft_fnc_A'*QA_init-...
            alfa_A/gama_A*xx_iA'.^gama_A*QA_init-...
            alfa_eA/gama_eA*xx_eA'.^gama_eA*QA_init+...
            pi_dom/(1-bet)*wq_A^(1-1/bet)*L_A*wage_dom_A*QA_init+...
            pi_frn_toUS/(1-bet)*wq_B^(1-1/bet)*L_A*flip(1-wage_dom_A)*QB_init+...
            wq_A*sum(QA_init)+Par.rev_US*trf_A*((1+trf_A)*(1+kapa)*eta/(1-bet))*wq_B*xk_B*flip(1-wage_dom_A)*QB_init;  
    IB_init = prft_fnc_B'*QB_init-...
            alfa_B/gama_B*xx_iB'.^gama_B*QB_init-...
            alfa_eB/gama_eB*xx_eB'.^gama_eB*QB_init+... 
            pi_dom/(1-bet)*wq_B^(1-1/bet)*L_B*wage_dom_B*QB_init+...
            pi_frn/(1-bet)*wq_A^(1-1/bet)*L_B*flip(1-wage_dom_B)*QA_init+...
            wq_B*sum(QB_init)+Par.rev_FN*trf_B*((1+trf_B)*(1+kapa)*eta/(1-bet))*wq_A*xk_A*flip(1-wage_dom_B)*QA_init; 

    IA_next = prft_fnc_A'*QA_next-...
            alfa_A/gama_A*xx_iA'.^gama_A*QA_next-...
            alfa_eA/gama_eA*xx_eA'.^gama_eA*QA_next+...
            pi_dom/(1-bet)*wq_A^(1-1/bet)*L_A*wage_dom_A*QA_next+...
            pi_frn_toUS/(1-bet)*wq_B^(1-1/bet)*L_A*flip(1-wage_dom_A)*QB_next+...
            wq_A*sum(QA_next)+Par.rev_US*trf_A*((1+trf_A)*(1+kapa)*eta/(1-bet))*wq_B*xk_B*flip(1-wage_dom_A)*QB_next;  
    IB_next = prft_fnc_B'*QB_next-...
            alfa_B/gama_B*xx_iB'.^gama_B*QB_next-...
            alfa_eB/gama_eB*xx_eB'.^gama_eB*QB_next+...
            pi_dom/(1-bet)*wq_B^(1-1/bet)*L_B*wage_dom_B*QB_next+...
            pi_frn/(1-bet)*wq_A^(1-1/bet)*L_B*flip(1-wage_dom_B)*QA_next+...
            wq_B*sum(QB_next)+Par.rev_FN*trf_B*((1+trf_B)*(1+kapa)*eta/(1-bet))*wq_A*xk_A*flip(1-wage_dom_B)*QA_next;  
    
    QwA_next = wage_dom_A*QA_next+flip(1-wage_dom_B)*QA_next*(1+kapa)/((1+kapa)*(1+trf_B))^(1/bet);
    QwB_next = wage_dom_B*QB_next+flip(1-wage_dom_A)*QB_next*(1+kapa)/((1+kapa)*(1+trf_A))^(1/bet); 
    
    gA_ss_vec(iter1) = (IA_next-IA_init)/IA_init;
    gB_ss_vec(iter1) = (IB_next-IB_init)/IB_init;
    
    if iter1>4/dt 
        
        g_rat = max(abs(((1+gA_ss_vec(end-4/dt))^(1/dt)-1)/...
                ((1+gA_ss_vec(end))^(1/dt)-1)-1),...
                abs(((1+gB_ss_vec(end-4/dt))^(1/dt)-1)/...
                ((1+gB_ss_vec(end))^(1/dt)-1)-1));
    
    end
    
    QQbar_rat(:,iter1) = [sum(QA_next)/QwA_next; sum(QB_next)/QwB_next];
    
    iter1 = iter1+1; 
    
end

gA_ss = mean(gA_ss_vec(end-2/dt:end));
gB_ss = mean(gB_ss_vec(end-2/dt:end));

qrat_Ass = QQbar_rat(1,end);
qrat_Bss = QQbar_rat(2,end);
Qrat_ss  = qrat_Ass/qrat_Bss;

madjX = +(round((1-bet)*(+1/bet*log((1+kapa)*(1+trf_B))-log(Qrat_ss))/log(lmd))-1);
madjX = min(max(madjX,-(mbar+1)),mbar); 

madjM = -(round((1-bet)*(-1/bet*log((1+kapa)*(1+trf_A))-log(Qrat_ss))/log(lmd))+1);
madjM = min(max(madjM,-(mbar+1)),mbar);

condX = abs(madj_X-madjX)>0;
condM = abs(madj_M-madjM)>0;

difff = 100*abs([psy*[gA_ss; gB_ss]/dt+ro-rr_guess(1:2);
                (QQbar_rat(:,end)-[qrat_A; qrat_B])]); 

disp(' ')
disp('   rr_guess    diff')
disp(' ')
disp([rr_guess difff(1:4)])
disp(' ')
disp('     m_X       m_M       qrat')
disp(' ')
disp([madjX -madjM Qrat_ss])
disp(' ')

varargout{1} = [vv_iA vv_iB xx_iA xx_iB xx_eA xx_eB];
varargout{2} = [prft_fnc_A/wq_A^(1-1/bet) prft_fnc_B/wq_B^(1-1/bet) wage_dom_A' wage_dom_B'];